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Abstract 

Smoothed particle hydrodynamics is a particle-based, fully Lagrangian, method for fluid-flow simulations. 
In this work, fundamental concepts of the method are first briefly recalled. Then, we present a thorough 
comparison of three different incompressibility treatments in SPH: the weakly compressible approach, where 
a suitably-chosen equation of state is used; and two truly incompressible methods, where the velocity field 
projection onto a divergence-free space is performed. A noteworthy aspect of the study is that, in each 
incompressibility treatment, the same boundary conditions are used (and further developed) which allows a 
direct comparison to be made. Problems associated with implementation are also discussed and an optimal 
choice of the computational parameters has been proposed and verified. Numerical results show that the 
present state-of-the-art truly incompressible method (based on a velocity correction) suffer from density 
accumulation errors. To address this issue, an algorithm, based on a correction for both particle velocities 
and positions, is presented. The usefulness of this density correction is examined and demonstrated in the 
last part of the paper. 
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1. Introduction 

Smoothed Particle Hydrodynamics (SPH) is a fully Lagrangian, particle-based approach for fluid-flow 
computations. This method was independently proposed by Gingold & Monaghan [T] and Lucy ,2J to 
simulate some astrophysical phenomena at the hydrodynamic level (compressible flow). Nowadays, SPH is 
more and more often used for flows with interfaces and in geophysical applications. Its main advantage over 
Eulerian techniques is that there is no requirement of the grid. Therefore, the SPH method is natural to 
use for complex geometries or multi-phase flows. 

An important issue in SPH is the proper implementation of the incompressibility constraints. In the 
present work, three different implementations are critically discussed and compared for a selection of val- 
idation examples. The first one is the Weakly Compressible SPH (WCSPH), which is the most common 
technique. It involves the standard set of governing equations closed by a suitably-chosen, artificial equation 



of state, cf. Sect. 4.1 The second and third implementations, called truly incompressible SPH (ISPH), are 
based on the Projection Method, introduced in SPH by Cummins & Rudman Generally in ISPH, a 
Poisson equation is solved to obtain the divergence-free velocity field, cf. Sect. |4.2| In the second approach, 
this Poisson equation is solved on an auxiliary regular mesh. Here, this method is called the Grid-Projected 
Poisson Solver (GPPS). Since this approach consists in combining the Eulerian and Lagrangian techniques, 
it is not an optimal choice for free-surface flows; yet, since SPH has a huge potential in the subject of 
multi-phase flows, its usefulness seems to be worth of discussion. The third approach considered is the 
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ISPH scheme with the Particle Poisson Solver (PPS), cf. Sect. 4.3 In this concept the Poisson equation is 
rewritten in SPH formulation and solved on the grid of Lagrangian points (particles) . The main advantage 
of this idea is that the additional domain discretization and related discretization errors are avoided. The 
paper is also related to the work of Lee et al. [3] , where a similar comparison between WCSPH and ISPH - 
PPS approaches has been done. However, in the present work, a comprehensive analysis of boundary condi- 
tions is performed. Moreover, since we use the same (ghost-particle) boundary conditions for both WCSPH 
and ISPH approaches, our comparison offers a new look at the usefulness of the incompressibility variants. 
For the sake of validation for both single- and two-phase flows, the three different implementations of the 
incompressibility constraints are compared for the lid-driven cavity flow and the Rayleigh- Taylor instability. 
For the former case, the impact of various computational parameters on the results have been analyzed and 
optimal setting has been found; it is consequently used also for the latter case. 

A very good comparison of ISPH solvers was presented by Xu et al. 0. However, the authors do not 
discuss the problem of accumulating density errors in the ISPH approach. Since, like other authors [3] [6], 
we have been faced with such a problem, we decided to examine the usefulness of the correction algorithm 
proposed by Pozorski and Wawrehczuk 7J, cf. Sect. [9] This formulation, used also in the PDF computations 
of turbulent flows [8 ■ , consists in retrieving the constant fluid density field by applying the correction to the 
particles' position. In practice, that procedure involves the computations of the second Poisson solver. The 
effectiveness of such an approach and the improvement of results are demonstrated. 

Another common field of research in SPH is the proper implementation of the boundaries. In the present 
work, among many techniques, the ghost-particle boundary approach is applied. This implementation 
involves the use of fictitious external particles that are reflections of the fluid particles in the computational 
domain. The main advantages of this method are simplicity and conformity with different phases of the 
fluid [H]. Computing the lid-driven cavity test problem with the WCSPH approach, we have found that the 
standard implementation of the no-slip boundary condition by the ghost-particle approach ^ may cause the 
stability problems. Our proposal to overcome these difficulties is based on the combination of the free-slip 
and no-slip boundary conditions. The detailed discussion of this problem is presented in Sect. [6j 

2. Formulation of the SPH method 

The main idea behind SPH is the introduction of the kernel interpolants for the field quantities so 
that the fluid dynamics is represented by particle evolution equations. The SPH method is based on three 
approximations . 

The first is interpolation of the field quantities at a point. To construct it, we utilize an integral 
interpolant A{r) of any field A(r) (for simplicity we consider here a scalar field) 

A{r)= [ A{r')W{r-r',h)dr', (1) 

where the integration is over all the domain D, and W{r, h) is a weighting function (the kernel) with the 
parameter h called the smoothing length (a linear dimension of smoothing). Generally, the kernel should 
posses a symmetrical form 

W{r,h) ^W{-r,h), (2) 



satisfy the limit condition 



lim W(r,/i) = (5(r), (3) 



where d{r) is the Dirac delta distribution, and should be normalized so that 

/ W{r,h)dr = 1. (4) 
Jn 

Additionally, the kernel should be at least as many times differentiable as the field A. Taking into consid- 
eration the computational effort and the proper implementation of the boundary conditions (Sect, pi), it is 



2 



worth to use compact support kernels. Since there are a lot of possibilities to choose W{r, h), we decided to 
compare solutions obtained employing three different 2D kernels (Sect. 7.1 1 : the cubic B-spline form 



1 - l?^ + h^: for < q < 1, 



10 



W{r,h) = ^{\{2~q)\ forl<g<2, (5) 



4 

0, otherwise, 



the quintic form proposed by Wendland (TU] 



^(r,,) 7 (1-1) (2, + l), for|r|<2., 
47rft.2 1^0, otherwise, 

and the quintic B-spline form introduced by Morris [llj 

{(3 -qf -6{2- qf + 15 (1 - qf , for < g < 1, 
{3-qf-6{2-q)\ forl<g<2, 
for2<,<3, 
0, otherwise, 

where q = \r\/h. 

The second approximation of the SPH technique is the discretization of space. It is done through dividing 
the domain into a fine-grained representation (particles). Each particle carries the properties of the field. 
The integral interpolant (•), Eq. ([l]), becomes then the summation interpolant (•) 

(A)(r)-^^(rfc)W^(r-rfc,/i)r!,, (8) 

b 

where rj, and J7b denote the position and volume of the particle b. The SPH task involves the foregoing 
computations of the interpolant at each particle, so that Eq. ([S]) may be rewritten into the common form 

{A)^^J2'^bWab{h)nb, (9) 

6 

where (A)^ = (A) (r,), A, = A(r,) and Wabih) = WUh) = W{n - Va, h). 

An additional advantage of SPH reveals with the differentiation of fields. In accordance with Q, the 
gradient of A{r) assumes the form 

V3(r) = / Vyl(r')VK(r-r',;i)dr'. (10) 



Taking advantage of the integration by parts rule and utilizing the kernel symmetry, we can transform the 
foregoing equation into 

VA{r)^ [ A{r')W{r~r')n'dS+ [ A{r')V'Wir-r',h)dr', (11) 
Jdn Jn 

where n' — n(r') is the normal vector to surface dfl. Generally, the first term does not necessarily vanish 
for finite domain sizes, cf. [12 . However, it is a common practice to neglect this term and deal with the 
boundaries explicitly. The SPH form (discretization) of (11) brings the common rule 

(^^)a = J2 AbyaWab{h)nb. (12) 
6 

Since the nabla operator acts only on the kernel, the gradient of the field is dependent only on the values of 
the fields at particles, not gradients. 



The third SPH feature is the assumption that the field value Aa at a point and its SPH approximation 
(A) are in relation 

(^>. « A,. (13) 

Since, in the case of the time dependent simulations, for each time step, the calculation of the field quantities 
is performed using the results obtained in the previous step, the above introduced assumption becomes 
natural. 



3. Governing equations 



The full set of governing equations for incompressible viscous flow is composed of the Navier-Stokes (N-S) 
equation 

^ = --Vp + j.V2u + f, (14) 
at g 

where g — const is the density, u the velocity, t the time, p the pressure, i> the kinematic viscosity and f an 
external force, and the continuity equation 



dg 
Itt 



= -gV • u. 



that for g = const arises to the form 



V • u = 0. 



(15) 



(16) 



The whole set of governing equations should be expressed in the SPH formalism. Utilizing the relation 

(17) 



(12 1, the divergence of velocity becomes 

(V-u)^ -^Ub- V,l¥a,(/i)l]b. 



Therefore, the continuity equation ( 15 ) takes the form 



dgg 
dt 



-Qa V Ufc • VaWab(h), 

. 9b 



(18) 



where m,})/ gi, — fi;,. It is important to note that various ways exist to express the divergence, for example 

dga _ 

b 



dt 



^mtUab ■ VaWafc(ft.), 



(19) 



where Uab = Uo — Uf,. The advantage of the above form over (18) is the symmetry with swapping particles 
a and b. Therefore, in practice, it is more accurate to use (19) |13j . However, there exists an alternative 
formulation. The fluid density can be computed directly from the SPH formula ([9| 



9a = ^gbWabWflb = ^TUbWabih). 



(20) 

b b 

A practical disadvantage of this approach is that g must be evaluated by summing over the particles before 
other quantities |TT]. Therefore, the computational effort increases (Sect. 7.1). Another disadvantage is 
the problem with representing sharp discontinuities near material interfaces. To avoid this difficulty, Hu & 
Adams [13] suggested to use the form 

9a=maJ2Wab{h). (21) 

b 
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In the SPH technique, the pressure N-S term is responsible for ensuring the incompressibiUty constraint 



(Sect.|4|). UtiUzing ([T2| it takes the form 



9 



Qa ^ Qh 



(22) 



Similarly to the continuity equation ( 19 1 there is a possibility to obtain more useful gradient approximations 

sed by 



In the present paper, we utilize the form proposed by Colagrossi and Landrini JT5 

Va +Pb 
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QaQb 



(23) 



is expressed as a combination of the finite difference and SPH discretisations [16] [17] . For the present work, 
we utilize the form (with a small regularising parameter rj — O.Ol/i) [17] 



(V(^V-u))^ = ^mJ8 



^ Vg + Vb Ugb ■ Tab 

Qa + Qb rli, + 772 



"^aWabih). 



Since SPH is fully Lagrangian approach, the particle advection equation completes the system 



(24) 



(25) 



4. Incompressibility treatment 

To compute incompressible flows in the Eulerian CFD, two approaches are commonly used: the ar- 
tificial compressibility method (where a specific equation of state is applied) and the pressure-correction 
technique (where the velocity field is projected onto the divergence- free space). Generally, the analogs of 
these techniques are also used in the SPH approach. They are presented in the following. 



4-1- Weakly Compressible SPH 

The most common technique is the weakly compressible SPH (WCSPH). It involves the set of governing 
equations closed by a suitably chosen, artificial equation of state p — p{g). Since the fluid pressure is an 
explicit function of g, the density gradient exerts an influence on the particle motion. The commonly used 
equation of state has the form [T^ 

o r / 

9 



P 



7 



1 



(26) 



where the reference density the numerical sound speed c and parameter 7 are suitably chosen to reduce 
the density fluctuations down to demanded level. In the present work, to assure density variations lower 
than 1%, we set 7 = 7 and c at the level at least 10 times higher than the maximal fluid velocity [TO] . 
However, since the sound speed is high, the time step should be, correspondingly, very small (cf. Sect. [5]). 
Therefore, the computational efficiency is the main weakness of the WCSPH method. 



4.2. Truly Incompressible SPH with the grid-projected Poisson solver 

The newest promising technique is the truly incompressible SPH (ISPH). It is based on the Projection 
Method which is a common approach for numerically solving time-dependent incompressible fluid-flow prob- 
lems. In this technique the pressure needed to ensure incompressibility is found by projecting the calculated 
velocity field onto the divergence-free space [3]. This is possible due to the Helmholtz decomposition which 
states: every vector field A that is twice continuously differentiable and vanishes faster than 1/r at infinity 
can be decomposed into the gradient and the curl parts as follows [20] 

A = V0 + VxB = Acurl free + Ajiv free, (27) 

5 



where </> and B are suitably chosen and 



V X Acurl free = V X (V0) = 0, 
V • Adiv tree = V • (V X B) = 0. 



(28) 



In the ISPH approach, the decomposition procedure begins with sphtting the integration of the N-S equation 
(14 1 on the time interval St = — into two parts. The first, so-called predictor step, gives the fractional 



velocity u* 



St 



2„ji 



(29) 



The right-hand side of the above equation contains all the N-S terms except the one connected with the 
pressure. The second part of the procedure is the correction step 



St 



— = --Vp 



n+l 



(30) 



It imposes the correction of u* to ensure compliance with the divergence- free constraint. To obtain an 
appropriate pressure we write the divergence of Eq. (301 



St 



(31) 



Since we expect a divergence-free velocity field at the end of the time step, we require that V • u"+^ — 0. 



Therefore, the formula (31) arises into the Poisson equation 

e 



St 



(32) 



Now, the correction step (30), performed withp"+^ obtained from the above relation, yields the divergence- 
free velocity field. 

The common way to solve the Poisson equation in the SPH approach is to compute it directly on 



particles (irregular grid), cf. Sect. 4.3 However, there exists another, less useful (specially for free-surface 



flows), but much more efficient treatment, called here the Grid-Projected Poisson Solver (GPPS). It consists 



in projecting the r.h.s. of (32) on an auxiliary grid and then to crunch it with some commonly known (from 
Eulerian approaches) solvers. This is a standard technique in various particle methods (more precisely, 
particle-mesh methods) such as Lagrangian PDF or Vortex-in-Cell. 

4-3. Truly incompressible SPH with the particle Poisson solver 

Employing together the SPH divergence and gradient operators, it is straightforward to obtain the 
direct SPH representation of the Laplace operator on the l.h.s. of (32 1. However, Cummins & Rudman 



[3] performed a simple one- dimensional hydrostatic test to show that such an approach produces a distinct 
pressure decoupling pattern. To avoid this problem, it is common to utilize the approximate Laplacian 



operator with the similar form as the viscous N-S term (24 1 



1. 



S-^rUb 4 Pab^ab ■ '^aWabih) 



b 



Qa Qa + Qb 



(33) 



where Pab =Pa—Pb- The pressure gradient is here computed using a finite-difference approximation. In this 



concept the Poisson equation ( 32 ) is solved on the irregular grid of Lagrangian points (particles) ; therefore, 
this variant of ISPH will be called here the Particle Poisson Solver (PPS). 

However, the main disadvantage of the PPS (Sect. |7.1[ ) is its inefficiency. Dealing with the Poisson 
equation using Eq. ( 33 ) consists in solving a linear equation system with a sparse irregular coefficient 



matrix. It requires much more CPU time and memory than the Poisson solvers performed on a regular grid. 
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5. Time step criteria 



In the way to assure the stabihty of the SPH scheme, several time step criteria must be satisfied [3] [TT] 
[T7] pT] . In the case of the ISPH approach, the CFL time step condition is 

^t<0.25 , ,^ , (34) 



where |u 



is maximal velocity in the flow. In the WCSPH approach, due to the utilization of the equation 



of state ( 26 1 , the CFL time step condition is 



St < 0.25 — . (35) 



Since we demand density fluctuations to be lower than 1%, we have chosen c > fO|u|maj; (Sect. 4.1). 
Therefore, when the flow is not dominated by viscous or external forces, WCSPH is computationally less 
efficient than the ISPH approaches. In the case of explicit schemes for viscous flows, another stability 
criterium is 

St < 0.125 — . (36) 
Additional condition must be satisfied due to the magnitude of particle accelerations f 

St < 0.25mm ( r^]' . (37) 

a \ \tn\ J 



6. Boundary conditions 

The proper implementation of the boundaries is one of the common topics in the SPH developments 
during recent years. Early stage applications of WCSPH involved high Reynolds number simulations with 
free-slip boundaries, performed using one layer of boundary particles placed at the wall. The layer exerted 
a strong repulsive force to prevent penetrating solid surfaces [22] . Since the number of interacting particles 
near the walls is decreased, the accuracy of numerical scheme degrades. Another treatment was proposed 



by Campbell [23] where the boundary condition was included in (111 through the residual boundary term. 
Today, the most often used boundary conditions are based on dummy particles [1] [23] . They are regularly 
distributed on the boundaries and have prescribed velocity during the whole simulation (no-slip condition). 



In ISPH, the Poisson equation ( 32 ) is solved on these dummy particles as well to repulse the fluid. To prevent 
inconsistency between density of inner particles and that of the wall, an additional set of dummy particle 
layers is placed outside the domain. Other popular, virtual-particle based boundary conditions utilize so- 
called mirror particles. These particles are given the prescribed velocity to assure the proper boundary 
condition. But, their properties are not integrated in time, unlike those of real SPH particles. Nowadays, 
there are two commonly used mirror-particle approaches. The first, developed by Morris [TT], consists in the 
combination of dummy and mirror particles. The velocity of inner particles is suitably projected on fixed 
boundary particles. Then, the boundary particles interact with the fluid. The second approach, so-called 
Multiple Boundary Tangent method, is similar to the previous one but the procedure of projecting particle 
velocities is different [25] . 

The mirror-particle techniques, presented above, have been further developed to another, more natural 
approach, i.e., the ghost particle method [3j. This technique is similar to the Classic Image Problem in 
electrostatics [20] [26]. To any particle a located at ra near the straight and infinite boundary, we introduce 
the image a' of this particle located at 

r,, = 2d + ra (38) 

where d is the vector pointing from the particle to the nearest point at the wall, cf. Fig.[T]^a). Since a chosen 
kernel is compact, the boundary may be finite. The role of these particles is to assure a high accuracy of the 
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computation (to remedy the lack of particles near the boundaries) and to enforce the boundary condition for 
the field quantities. Thus, the natural way to obtain the proper implementation of the boundary condition, 
with the no-slip condition for velocity field, is to set: 

Ua'=2Uw-Ua, Vfla' ^ nia, Qa' = Qa, (39) 

where is the velocity of the boundary. To enforce the proper Neumann boundary condition for the 
pressure 



where n is the vector normal to the boundary, we extend the set ( 39 ) by an approximation of ( 40 ) 

Pa'^Pa- (41) 

Let stand for a position of any point placed at the boundary. Due to the kernel symmetry we have 

Wir^-ra,h)^W{r^~ra',h). (42) 
Now, utilizing the SPH summation interpolant ([9| for velocity, we write 

—UaW{r^ - Ta, h) + } — ^(2u„ - Ua)W{r^ - r^, h) = 

Ba ^ Qa' 

m„ 



2uw —W{v^ ~ra,h) 

^ Ba 



(43) 



a 



Since the summation in the above equation is over fluid particles only, assuming nearly homogeneous dis- 
tribution of fluid particles, we have 

Y^^W{r^-v^^h)^\ (44) 



Ba 



and (43) becomes 

(u) (rw) w Uw. (45) 



Therefore, (39) provides the proper formulation of the no-slip condition. 

Another boundary type that can be treated with the ghost-particle approach is an inner corner. The 
technique of constructing particle images is presented in Fig. [ijb). In this case we have to use three mirror 
particles. It is important to note that the influence range of this corner is smaller than 2h. For a larger 



2h 



2h 



2h 



a=-l 



0=0 



a=l 



2h 



2h 



Cummins like 



proper boundary 
condition 



Figure 2: Tlie problem with the ghost-particle no-shp boundary condition near the corner: depending on value of the parameter 
a, the velocity boundary condition changes. 



distance from the corner, the boundary condition boils down to the previous case (both in vertical and 
horizontal directions). For the distances smaller than 2h, Cummins and Rudman [3] use the third particle 
placed symmetrically in respect to the corner point possessing the same density and mass but opposite 
velocity to fluid particle. Since the velocity vector computed at the corner is non-zero in the case of not 
moving boundaries, this formulation is improper. Therefore, in the way to find a more accurate approach, 
we parametrize the mirror particles' properties as follows (using the notation from Fig. [T]): 



Qa Qa' Qa" 

ma = TO„' = m„ 



Qa'' 



Via" = -Ma, Va" = 2d' + Ta, 

Ua'" = 2aUw + Ua, Ta'" = 2d + 2d' + Tq, 



(46) 



where a may change from —1 up to 1. Now, for a point at the boundary, we may carry similar investigation 



as in ( 43 ) 



(u) (rw) = V —ViaW{v^ - Fq, ft,) + V — (2Uw - \la')W{v^ - Ya'.h) 

+ V —M,"W{v^ - v^u^h) + V ^(2uwa + vl':)W{v^ - r,,., ft) 



(47) 



= 2u„ ^ — [W{v^ - 2d ra, ft) + aW{v^ - 2d - 2d' + r,,, ft)] . 

a 

For a point placed exactly at the corner, the above equation reduces to the form 

1 



(u) (re 



-Uw(l + a). 



(48) 



Depending on value of the parameter a, the velocity of the corner changes from u(rcorncr) = (for a = —1) 
up to u(rcorncr) — Uw (for a — 1) . Despite this parametrization, there is no way to assure proper velocity 
at boundary that is closer than 2ft from the corner. This issue is presented in Fig. [2] Yet, the simulations 
of various test problems show that varying parameter a does not have a significant impact on the global 
(except near the corner) velocity or density fields. Therefore, in the present work, we utilize the no-slip 
condition with a = 0. 
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X X 



Figure 3: The WCSPH results of the hd-driven cavity {Re = 1000) at t = 0.03 with: a) the no-shp and b) the free-slip boundary 
treatment for velocity divergence computation; employing no-slip condition induces instabilities near the corners. 



Another problem with the no-sUp condition implementation appears during computation of density in 



WCSPH utilizing the continuity equation ( 15 1 . Assuming a statistically homogeneous distribution of particle 



positions, the velocity component tangential to the boundary is negligible for computing the divergence. 
Therefore, along straight boundaries the no-slip condition is properly stated. The problem appears near the 



corners. Utilizing (12 1, we may compute the divergence of velocity at the boundary near the corner 



(V • u) (r^) = — • Va^(rw - r,, h)+J2 — Ua' • ^a'Wir^ - v,,,h) 

On Qa' 



Qa 

(49) 

V —Mia" ■ Va.Ty(rw - Va",h) + V "^Ug,,, ■ Va'"W{v^ - Va"',h). 
^ Qa" ^ 0„n, 



For the point Tc that is placed in the fluid e — > away from the corner we have 

VaW{Yc-Va,h) ■ Xi^ ^ V a'W {v ^ - V a' ,h) -11^ = 
-Va"VF(rc - Va",h) ■ = -Va"'W{Vc - Va"',h) ■ n^r 



and 



VaW(rc - Ta, h) ■ riy = -Vq' W(rc - Ya-.h) ■ riy = 

Va"W{Vc - Va",h) ■ THy = -Va"'W{Vc - Va"',h) ■ Uy, 



(50) 



(51) 



where n^. and n^are unit vectors in x and y directions respectively. Therefore, connecting the above relations 
with ( 49 ) and ( 46 1 we obtain 

(V • u) (r,) « 0. (52) 

Let us consider the situation where fluid particles are driven to the corner. To prevent penetrating the 
boundary, there should appear a repulsive force near this corner. In WCSPH this is done by a local density 
increase. However, since the divergence of velocity is always close to zero near the corner, according to the 



continuity equation (15), the density does not change. This may induce a growth of velocity instabilities 



near the corners (see Fig.[3|). Therefore, only for computing the divergence of velocity, we suggest to use the 
free-slip condition. For a straight boundary, cf. Fig.|4j the free-slip mirror particle is placed at r^' ~ 2d + ra. 
It has the same mass and density, its velocity component normal to the boundary is opposite, while the 
tangential component is unchanged. At the corner, cf. Fig. |4];b), the mirror particles carry the following 
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wall 



fluid 



Figure 4: The ghost-particle free-slip boundary treatment for: a) the straight wall and b) the corner. 




properties: 



Qa 

ma 



ga' 

rria' 
2d 



Ta" = 2d' 
Ta'" = 2d - 



= Qa" 
= ma 

2d' 



Qa" 



u 



a 



a" 



nx 



-Ua ■ Ux 
~Ua ■ Ux 



U„ 



-Ua • n 
Ua ■ n 



J" 



^a 



•■V 



(53) 



Obviously, this kind of boundary treatment presented here is basicaUy hmited to flat wall segments only; 
more advanced formulations exist and represent the state of the art for arbitrarily-shaped boundaries, e.g., 
|12j . However, since the main focus of our paper is on incompressibility treatments, we have chosen simple 
geometry cases to have enough reference data. 



7. Velocity error measurement 

7.1. Lid-driven cavity 

The lid-driven cavity is a common test of numerical algorithms for viscous flows. It involves a fluid at 
density Qq inside a square (L x L) box where only one boundary moves with the constant velocity n^. The 
geometry is very simple, however there is no analytical solution. In the present work, we computed the 
lid-driven cavity flow at Re = |u^|L/j/ = 1000. For this value of Reynolds number the flow is still laminar 
and there is no necessity to use a turbulence model. All results are suitably non-dimensionalised with L, 
|uu,|, qq (especially time is normalized with the convective time scale L/|utu|) and compared to those from 
a numerical calculation on a fine grid performed with the Eulerian solver by Ghia et al. |27j . 



7.1.1. The kernel type influence 

One of the important issues that has an impact on the SPH solutions is a proper selection of the kernel. To 
compare the treatments of the incompressibility constraint, we decided to perform the benchmark simulation 
of the lid-driven cavity flow with the three kernels presented in Sect. [2] The particles' number has been 
chosen with two requirements: N should be large enough to get solutions comparable to reference data (yet 
not fully converged); however, it should be sufficiently small to show kernels' defects. We decided to use 
N = 3600 particles in the domain and h/Ar — 2, cf. Sect. 7.1.2 The steady-state solution velocity profiles 
are presented in Fig. [5] The results performed with the cubic spline kernel are the most inconsistent with the 
reference data. This discrepancy is caused by the particles' clustering phenomenon. Figs.[6]and[7]respectively 
present particles' spatial distribution and histograms of the distance between the nearest pairs of particles 
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for all mentioned kernels. As it transpires from the histogram for the cubic spline kernel, in this case there 
are many particles that move joined in groups, cf. also details in Fig. [6j as a consequence, the accuracy of 
the scheme is radically decreased. As far as the quintic kernels are concerned, both perform similarly to each 
other. However, the kernel proposed by Morris et al. ^ shows a more pronounced tendency to clustering 
(Fig. [7]); moreover, it is not zero up to \r\/h = 3 (as contrasted to \r\/h = 2 for the other kernels), so it 
involves many more particles. Therefore, the computational time is increased about 12% comparing to the 
Wendland kernel. Summarizing, due to a good agreement with the reference data, not noticeable clustering 
phenomenon, and finally the efficiency, for the further analysis we decided to use the quintic Wendland 
kernel. An interesting work about its behavior in SPH has recently been performed by Robinson [28j (Ch. 
7). 

7.1.2. The kernel size influence 

Apart from the initial inter-particle distance Ar (related to the number of particles in the domain), in 
the SPH approach there exists another characteristic length - the kernel parameter h. Examining the impact 
on the solutions, we decided to compute the lid-driven cavity problem {Re — 1000) for chosen N = 3600, 
the Wendland kernel and different h/Ar values. Figure |8] presents the velocity profiles for h/Ar = 2.31, 
2.0, 1.67 and 1.5 computed with all the SPH incompressibility variants. In comparison to the Ghia et 
al. reference data, there is no significant effect of the kernel size on the velocity field. However, since the 
parameter h/Ar determines the number of particles under the kernel hat, the number of interaction between 
particles (and the computational effort) grow like ~ (h/Ar)^, where D is the space dimension. Weighting 
between the computational times and the accuracy, for all simulations presented henceforth, we decided to 
use h/Ar = 2.0. 

7.1.3. The particle number influence 

Examining the influence of the spatial resolution, the simulations were performed with a different number 
of particles in the domain (from N — 1600 up to iV = 57600). The steady-state velocity profiles are presented 
in Fig. [9j In all the methods of incompressibility treatment, the profiles obtained for N = 1600 did not 
agree very well with the Ghia solution [27] obtained on 129 x 129 Eulerian mesh. Performing simulations 
with the higher resolutions, the accuracy increases so that in the case of iV = 14400, independently on the 
compressibility treatment, the obtained results are in very good agreement with the reference simulations. 
For N — 57600 the velocity profiles computed using WCSPH and ISPH - PPS approaches practically overlap 
with the Ghia's reference data. Unfortunately, the ISPH - GPPS solution is not so accurate as we expect. 
This deficiency is due to a higher numerical diffusion caused by the projection of the quantities from particles 
on a regular grid. 

Due to the utilization of the equation of state and the time step constraint in the WCSPH approach, 
the CPU time of the lid-driven cavity simulation computed by the ISPH method with PPS is about 7 
times shorter. On the other hand, performing the ISPH-GPPS simulations, the computational effort may 
be reduced about 15 times. The comparison of CPU times for all considered SPH schemes is presented in 
Fig. [101 

7.2. The Ray leigh- Taylor instability 

The Rayleigh- Taylor instability is one of generic multi-phase flows, therefore it is commonly utilized as a 
test problem. It involves two immiscible fluids enclosed in a rectangular domain of width L and height 2L. 
Initially, in our case, the phases are separated by the interface located at y = 1 — 0.15 sin {2Trx). The lower 
component has density ql = go, while the upper one gjj = l-Sgg. Since the system is subject to gravity 
g — (0, ~g) and the upper phase is heavier, in the absence of the surface tension an instability always arises 
and vorticity is generated. The Reynolds number may be defined as 

Re = = 420. (54) 

v 

The simulations were performed for 120 x 240 particles in the domain. To compute the hydrostatic force, we 
use the technique described in 29J , where the hydrostatic pressure is computed on a regular mesh and later 
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Figure 5: The lid-driven cavity steady-state velocity profiles for: (a) WCSPH, (b) ISPH-GPPS and (c) ISPH-PPS against Ghia 
et al. |27| results; profiles obtained for different kernels; A' = 3600, h/ Ar = 2. 
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Figure 6: The lid-driven cavity steady-state solution {Re = 1000, N = 3600) computed with the WCSPH approach and kernels: 
fsjl, the particle clustering phenomenon is noticeable with the cubic spline kernel. 
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Figure 7: Histograms of the distance between the nearest pairs of particles; the results obtained for the WCSPH approach and 
kernels: (jej, |7|. 
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Figure 8: The lid-driven cavity steady-state velocity profiles for: (a) WCSPH, (b) ISPH-GPPS and (c) ISPH-PPS against Ghia 
et al. |27| results; profiles obtained with diflFerent /i/Ar values; results obtained using the Wendland kernel |10) and A'^ = 3600 
particles in domain. 
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Figure 9: The lid-driven cavity steady-state velocity profiles for: (a) WCSPH, (b) ISPH-GPPS and (c) ISPH-PPS against Ghia 
et al. |27| results; results for different number of particles N; data obtained using the Wendland kernel 10 and h/Ar = 2. 
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Figure 10: The CPU times to obtain the steady-state solution of the Ud-driven cavity {Re = 1 000) using the WCSPH and both 
ISPH approaches; for the WCSPH method, two continuity equations are compared: Eq. l |19[ l and density definition \21\ . 



projected on the particles. Figure 11 presents particle positions (directly showing the location of liquid- 
liquid interface) computed with different treatments of the incompressibility condition. All the simulations 
are compared to the reference solutions from the Level-Set method (312 x 624 cells) computed by Grenier 
et al. [30]. Presented data were obtained at t = b (normalized with the convective time scale y/L/g). 

Comparing the interface shapes, we find no convincing arguments to judge which SPH approach is more 
accurate. However, the comparison of the computational times shows that the most valuable choice is the 
ISPH-GPPS. The ISPH-PPS is slower about 6 times, while the WCSPH about 13 times. 



8. Density error measurement 

In order to measure the density estimation errors, we compute two quantities: the mean density over 
the flow domain 'g{t) and the root mean square of the density fluctuations Qrmsit)- However, since in the 
Projection Methods the values of the density carried by particles do not change, there is the necessity to 
use formulas which are able to compute 'g(t) and Prms(i) taking into account both the density values and 
spatial distribution of the particles. The simplest proposals (expressed in the SPH formulae) are: for the 
mean density 



1 



ma 

b 



(55) 



and for the density fluetuations ^'rms(^) 



a \ b 



m 



(56) 



In WCSPH, comparing the influence of the kernel shape, we observe that the quintic kernel proposed 
by Wendland ([6| gives the smallest fluctuations of the density field, cf. Fig. [T2{^a). The weakness of the 
Wendland kernel is a slight overestimation of the mean density. Even for a homogeneously distributed set 
of particles {t = 0) the density field is flawed, cf. discussion in [31 . At decrease of the h/Ar parameter 
causes this error to grow from 0.2% for h/Ar — 2.31 up to about 0.8% for h/Ar — 1.5, cf. Fig. [l2jb). 
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Figure 11: The Rayleigh- Taylor instability; particle positions at t = 5 obtained using different incompressibility treatments; 
solid line: liquid-liquid interface from the reference Level-Set solution |30j . 



Interestingly, parameter /i/Ar has no significant effect on g,. 
with WCSPH is presented in Fig. 12 'c). The density r.m.s. 
domain, stabilizes after about t = I and remains less than 1% of initial density. Moreover, with increasing 
iV the density r.m.s. goes down. In the case of ISPH solvers, the growth of the r.m.s. stops only after 



is(t). The particle number influence obtained 
independently on the number of particles in 



obtaining the steady-state solution (roughly at t = 120), cf. Fig. 12 'd). Intriguing is the fact that for the 



ISPH approach, the increase of the number of particles N increases the r.m.s. For ISPH-PPS the r.m.s. of 
the steady-state solution changes: from 7% for N = 1600 up to 9% for N = 14400. For ISPH-GPPS the 
density field is additionally affected by the projections between the particles and the grid. These projections 
cause additional smoothing and the density field r.m.s. of the steady-state solution varies: for N — 14400 
about 7%, for N = 3600 about 6% and N = 1600 about 5% of initial density. 

This discrepancy between WCSPH and ISPH density error values is caused by the completely different 
mechanisms of computing the dynamical pressure in the Navier-Stokes equation. Since in the case of the 
WCSPH approach, the suitably-chosen equation of state is used, the density field suffers from short-scale 
density waves presented in Fig. [l3j[a). In the case of ISPH, there is no explicit correction to the density. 
Therefore, during the computations, the density field strongly deviates from the level of the initial state, cf. 



Fig. 13 'b). Since the Projection Methods give such considerable density errors, it is meaningful to implement 



and an additional correction procedure, cf. Sect. [9j 

In the case of the Rayleigh- Taylor instability, the comparison of density errors (cf. Fig. 14 ) shows similar 
behavior as in the lid-driven cavity case. 



9. Constant-density ISPH approach 

In the previous section we have shown that in the ISPH approach the density errors cumulate during 
computations. To solve this problem, Pozorski and Wawrehczuk ^ suggested to use the second correction 
term (a similar approach was later proposed, in the multi-phase flow context, by Hu and Adams Let 
us consider incompressible fluid with initially uniform density g^. After the time step performed with the 
ISPH approach, the corrected velocity field satisfies the divergence-free condition; however, this procedure 
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Figure 12: The density mean value and the r.m.s. obtained for the lid-driven cavity Row at Re = 1000; the effect of: (a) the 
kernel choice, (b) h/Ar, (c) number of particles A'^ influence in the WCSPH approach; (d) particles number N influence in 
both: ISPH-PPS and ISPH-GPPS techniques. 
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Figure 14: The density mean value and the r.m.s. obtained for the Rayleigh- Taylor instability (only upper phase qjj = l.Sgo) 
using WCSPH and both ISPH approaches. 
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Figure 15: The density field with local disturbance; left: the initial state (regular set of particles with one particle displaced), 
right: after 50 correction iterations. 



does not explicitly guarantee that g — const. The Pozorski and Wawrenczuk idea consists in performing an 
additional correction to the particle positions 



1 



u: 



n+l 



5t 



1 



(57) 



where u"+^ and f""''"'^ are respectively the divergence-free velocity field and particle positions obtained after 



the ISPH time step, Eqs. (30) and (32 1, while p* appears as a potential correction field computed from 



(58) 



where g^'^^ is the density obtained after the ISPH time step. Both and g^'^^ are computed using Eq. (21 ) 



The Poisson Eq. ( 58 ) is obtained from the request 



b 



-\h) 



b 



i.n+1 
'■b 



^P*b.h) = g^iva 



(59) 



by the Taylor series expansion around — f^^^. 

Since, performing such a correction, the second Poisson equation has to be solved, the computational 
effort is increased. However, there is no necessity to compute the correction at each time step; rather, it is 
applied only if the density error exceeds a certain threshold value. On the other hand, the procedure may 
be performed several times in one time step. Figure [15] shows how the initial disturbance of the density field 



(a regular set of particles with one particle displaced) is corrected with ( 58 ) 



To perform the test problems, we decided to utilize ISPH with particle positions correction performed 
once per time step. Since both Poisson equations were solved on particles, we got fully grid-free approach. 
The velocity profiles of the lid-driven cavity problem [Re 



1000) are presented in Fig. 16 



There are no 

noticeable differences between the velocities obtained using ISPH solver with (Fig. 16 1 and without (Fig. |9| 
the Pozorski and Wawrenczuk correction. However, comparing the density fields, cf. Fig. [17) with the results 



obtained without the correction, cf. Fig. 13 d), the advantages of this approach become obvious: the growth 
of the r.m.s. density stops after t = 5 and stays at a level less than 0.1%. The regularisation of the particles' 
distribution is also visible in the histograms of the distance to the nearest neighbor of each particle: the 



comparison of histograms for ISPH-PPS with and without the correction procedure is presented in Fig. 18 
Another convincing argument for the utility of the Pozorski and Wawrenczuk approach is the computational 
effort. Since, as in the case of WCSPH, the density error does not accumulate, and, despite the use of two 
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Figure 16: The lid-driven cavity velocity profiles at the steady-state against Ghia et al. [27 ] reference data (fie = 1000); results 
obtained using ISPH with PPS and Pozorski & Wawrenczuk correction [7] for the Wendland kernel 10 , h/Ar = 2 and different 
number of particles N. 



Poisson solvers, the computational time is about 3 times shorter (and about 2 times longer than ISPH-PPS), 
the use of the density correction algorithm seems to be profitable. 

On the other hand, it is interesting to note a higher convergence rate with increasing number of particles 
in domain, as compared to the solutions obtained by Lee et al. [3| (Fig. 6) and Xu et al. [S] (Fig. 25). We 
suppose that a better convergence rate in our case is due to a proper choice of computational parameters 
[h/Ar, kernel type and b.c). 

For the Rayleigh- Taylor instability, the calculated particle positions at times i = 1, 3 and 5 are presented 



in Fig. 19 The comparison with the reference interface shapes fSff show quite a good agreement. What is 
more, there are no considerable differences between ISPH with and without the Pozorski and Wawrenczuk 
correction term. Moreover, as in the case of the lid-driven cavity, the use of such a correction term assures 
that the density error (Fig. 17) is kept below the desired level and the computational time is still acceptable 
(in comparison with the WCSPH approach). 



10. Conclusions 

In this paper three different SPH incompressibility treatments were considered. To validate these ap- 
proaches, two tests cases were simulated: the lid-driven cavity at Re = 1000 and the Rayleigh- Taylor 
instability at Re = 420. Summarising, in comparison to the previous tests presented by Lee et al. [J and 
Xu et al. [5], in all of the SPH incompressibility treatments, the velocity field exhibit a better convergence 
rate, when the number of particles in the domain grows. Detailed study of the incompressibility treatments 
revealed some pros and cons of the approaches. The main disadvantage of the WCSPH approach is its high 
computational cost. Therefore, from this point of view, it is better to use the ISPH approaches. Neverthe- 
less, both the ISPH-GPPS and the ISPH-PPS solvers in their original formulation suffer from the density 
error accumulation. 

For the fluid flows where besides the properly modelled velocity field, the correct density field is of 
importance, an interesting alternative to WCSPH is the ISPH approach with the separate density correction. 
In this approach the density accumulation error can be reduced to a specified level while still retaining the 
grid-free formulation. It is important to note that, as such, none of the ISPH approaches can be used to 
properly model gas-liquid flows such as a bubble raising in water, where the divergence-free condition is 
satisfied only for the liquid phase. Such flows can be relatively easily modelled with the WCSPH approach, 
where each phase can be computed with a different equation of state. 
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Figure 17: The density mean value and r.m.s. obtained using the ISPH approach with PPS and Pozorslii &; Wawrehczuli density 
correction [7]; (a) the lid driven cavity [Re = 1000), (b) the Rayleigh-Taylor instability [Re = 420); the results obtained for 
the Wendland kernel JO] and h/Ar = 2. 
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Figure 18; Histograms of the distance between the nearest pairs of the particles; the results obtained for the ISPH-PPS approach 
(a) without and (b) with the Pozorski and Wawrehczuk correction procedure. 
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Figure 19: The Rayleigh- Taylor instability (Re = 420) computed using ISPH with PPS and Pozorski & Wawrenczuk density 
correction [3; the black lines denote the interface position obtained by Grenier et al. |30j using the Level-Set formulation 
(312 X 624 cells). 

The study has also shown the importance of correct choice of computational parameters: the kernel 
formula W, its smoothing length h for given inter-particle distance h/Ar, and the number of particles N 
in the domain. The impact of those parameters on the results has been analyzed and optimal kernel has 
been found (the Wendland formula). Regarding the other quantities, they represent a compromise between 
the CPU efficiency and accuracy. Further work will be undertaken to extend the ISPH approach making it 
possible to model gas-liquid mixed flows with the interface phenomena such as the surface tension or the 
Marangoni effect. Also, there is the need of appropriate wall b.c. for SPH simulations in complex geometries 
(arbitrary shape of the boundary). 
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List of figure captions 



Fig. 1: The ghost-particle no-shp boundary scheme: a) straight wall, b) inner corner. 

Fig. 2: The problem with the ghost-particle no-slip boundary condition near the corner: depending on 
value of the parameter a, the velocity boundary condition changes. 

Fig. 3: The WCSPH resuhs of the lid-driven cavity {Re = 1000) at t = 0.03 with: a) the no-slip and b) 
the free-slip boundary treatment for velocity divergence computation; employing no-slip condition induces 
instabilities near the corners. 

Fig. 4: The ghost-particle free-slip boundary treatment for: a) the straight wall and b) the corner. 

Fig. 5: The lid-driven cavity steady-state velocity profiles for: (a) WCSPH, (b) ISPH-GPPS and (c) 
ISPH-PPS against Ghia et al. [22] results; profiles obtained for different kernels; N = 3600, h/Ar = 2. 

Fig. 6: The Ud-driven cavity steady-state solution {Re = 1000, N = 3600) computed with the WCSPH 
approach and kernels: ([t]); the particle clustering phenomenon is noticeable with the cubic spline 

kernel. 

Fig. 7: Histograms of the distance between the nearest pairs of particles; the results obtained for the 
WCSPH approach and kernels: 0. 

Fig. 8: The hd-driven cavity steady-state velocity profiles for: (a) WCSPH, (b) ISPH-GPPS and (c) 
ISPH-PPS against Ghia et al. |27j results; profiles obtained with different h/Ar values; results obtained 
using the Wendland kernel ITO] and N = 3600 particles in domain. 

Fig. 9: The hd-driven cavity steady-state velocity profiles for: (a) WCSPH, (b) ISPH-GPPS and (c) 
ISPH-PPS against Ghia et al. [27] results; results for different number of particles N; data obtained using 
the Wendland kernel [10] and h/Ar — 2. 

Fig. 10: The CPU times to obtain the steady-state solution of the lid-driven cavity {Re — 1000) using 
the WCSPH and both ISPH approaches; for the WCSPH method, two continuity equations are compared: 
Eq. (19) and density definition (21). 

Fig. 11: The Rayleigh- Taylor instability; particle positions at i = 5 obtained using different incompress- 
ibility treatments; solid line: liquid-liquid interface from the reference Level-Set solution |30j . 

Fig. 12: The density mean value and the r.m.s. obtained for the lid-driven cavity flow at Re = 1000; the 
effect of: (a) the kernel choice, (b) h/Ar, (c) number of particles N influence in the WCSPH approach; (d) 
particles number N influence in both: ISPH-PPS and ISPH-GPPS techniciues. 

Fig. 13: The lid-driven cavity flow: density field at the steady-state solution. 

Fig. 14: The density mean value and the r.m.s. obtained for the Rayleigh- Taylor instability (only upper 
phase gu = l.Sgo) using WCSPH and both ISPH approaches. 

Fig. 15: The density field with local disturbance; left: the initial state (regular set of particles with one 
particle displaced), right: after 50 correction iterations. 

Fig. 16: The lid-driven cavity velocity profiles at the steady-state against Ghia et al. [57] reference data 
{Re = 1000); results obtained using ISPH with PPS and Pozorski & Wawrehczuk correction [7 for the 
Wendland kernel [10], h/Ar = 2 and different number of particles N. 

Fig. 17: The density mean value and r.m.s. obtained using the ISPH approach with PPS and Pozorski & 
Wawrehczuk density correction [7]; (a) the lid driven cavity {Re = 1000), (b) the Rayleigh- Taylor instability 
{Re = 420); the results obtained for the Wendland kernel [T^ and h/Ar = 2. 

Fig. 18: Histograms of the distance between the nearest pairs of the particles; the results obtained for 
the ISPH-PPS approach (a) without and (b) with the Pozorski and Wawrehczuk correction procedure. 

Fig. 19: The Rayleigh- Taylor instability {Re = 420) computed using ISPH with PPS and Pozorski & 
Wawrehczuk density correction [7]; the black lines denote the interface position obtained by Grenier et 
al. [5D] using the Level-Set formulation (312 x 624 cells). 
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